/*
 * KineticProblem.cpp
 *
 *  Created on: 18 Oct 2011
 *      Author: allan
 */

#include "KineticProblem.h"
#include "Utils/Assert.h"

KineticProblem::KineticProblem()
{}

KineticProblem::KineticProblem(const Multiphase& multiphase, 
                               const ReactionSystem& eReactions, 
                               const ReactionSystem& kReactions, 
                               const vector<string>& jSpecies, 
                               const vector<string>& eSpecies, 
                               const vector<string>& kSpecies)
{
	Initialiaze(multiphase, eReactions, kReactions, jSpecies, eSpecies, kSpecies);
}

void KineticProblem::Initialiaze(const Multiphase& multiphase,
                                 const ReactionSystem& eReactions,
                                 const ReactionSystem& kReactions,
                                 const vector<string>& jSpecies,
                                 const vector<string>& eSpecies, 
                                 const vector<string>& kSpecies)
{
	// Set data members
	this->multiphase = multiphase;
	this->eReactions = eReactions;
	this->kReactions = kReactions;
	
	// Set the number of primary (Nj) and kinetic (Nk) species
	Nj = jSpecies.size();
	Nk = kSpecies.size();
	
	// Set the indexes of the kinetic species
	kIndexes = multiphase[kSpecies];
	
	// Assemble the stoichiometric matrices (ukj) and (ukk)
	ukj = kReactions.AssembleStoichiometricMatrix(jSpecies);
	ukk = kReactions.AssembleStoichiometricMatrix(kSpecies);
	
	// Initialiaze the total equilibrium solver (eSolver)
	eSolver.Initialiaze(multiphase, eReactions, jSpecies, eSpecies);
}

void KineticProblem::SetParameters(double T, double P, const VectorXd& n)
{
	this->T = T;
	this->P = P;
	this->n = n;
}

void KineticProblem::Function(double t, const VectorXd& ujnk, VectorXd& R)
{
	// Transfer the uj and nk part of the vector (ujnk) = [uj nk] to the vectors (uj) and (nk) respectively
	const VectorXd uj = ujnk.segment(00, Nj);
	const VectorXd nk = ujnk.segment(Nj, Nk);
	
	// Update the nk part of the vector (n)
	for(unsigned k = 0; k < Nk; ++k) n[kIndexes[k]] = nk[k];
	
	// Solve the equilibrium problem with total composition of the primary species (uj)
	auto diagnostic = eSolver.Solve(uj, T, P, n);
	
	Assert(diagnostic.converged, "Method did not converged.");
	
	// Calculate the activities of the species, now that the composition vector (n) has been determined
	const VectorXd a = multiphase.Activities(T, P, n);
	
	// Calculate the mass rates of the kinetic-controlled reactions with temperature (T), pressure (P), composition vector (n), and activity vector (a)
	const VectorXd rk = kReactions.KineticRates(T, P, n, a);
	
	// Calculate the right-hand side of the ordinary differential equations
	R.segment(00, Nj) = ukj.transpose() * rk;
	R.segment(Nj, Nk) = ukk.transpose() * rk;
}
